module define
     parameter (N_SGL_DIGITS=6)
     parameter (N_DBL_DIGITS=15)
     parameter (N_INT_ORDER=8)
     parameter (NBS=selected_real_kind(N_SGL_DIGITS))
     parameter (NBD=selected_real_kind(N_DBL_DIGITS))
     parameter (NBI=selected_int_kind(N_INT_ORDER))
!#ifdef DBL         
     parameter (NBF=NBD)
!#else         
!         parameter (NBF=NBS)
!#endif         
!
! constants
    parameter (COL_AB = 1)
    parameter (COL_SC = 2)
    integer ng2,ndecgrp
    logical TRUE, FALSE                                         
    real(NBF)  cinf,zero,one,half,big,CKELVIN,PI
    parameter (ng2=2,ndecgrp=6)
    parameter (TRUE=.true.,FALSE=.false.)

    parameter (CKELVIN=273.15_NBF,PI=3.1415927_NBF)
    parameter (cinf=1.0e30_NBF,zero=0._NBF,one=1._NBF,EPSILON=1.e-30_NBF,half=0.5_NBF,big=1.0e30_NBF)

    integer(NBI) UP,DOWN,CIRCLE
    parameter (UP=1,DOWN=2,CIRCLE=3)

    integer, parameter :: HMZMC = 1, HMZRT = 2
    type scatmat
        sequence
        integer(NBI) ib,ie
        real(NBF),pointer :: from(:)
        real(NBF),pointer :: to(:)
    end type
    type cmtofm
        sequence
        integer(NBI) nfm
        integer(NBI),pointer :: fm(:)
    end type

    character(132) mesg
    character(20) version

!   for mc simulation
    type particle
        sequence
!        real(8) :: xy(2), dir(2), mu, w
!        integer(4) :: lm(2), g, dummy
!        integer(4) :: seed(2)
        real(8) :: xy(2), w
        integer(4) :: lm(2)
    end type
    integer, parameter :: DW=1, DE=2, DN=3, DS=4
    integer, parameter :: CENT=1, EDGE=2
    
    ! for larsen's CMFD =====================================================
    integer, dimension(4) :: LTSIDXX ! surface index x-dir
    data LTSIDXX /-1,+1,-1,-1/
    integer, dimension(4) :: LTSIDXY ! surface index y-dir
    data LTSIDXY /-1,-1,+1,-1/
    integer, dimension(4) :: LTPHI2X ! make Omega_x positive
    data LTPHI2X /-1,+1,0,0/
    integer, dimension(4) :: LTPHI2Y ! make Omega_y positive
    data LTPHI2Y /0,0,+1,-1/
    integer, dimension(4) :: LTSURFDIR ! 1 for x-dir, 2 for y-dir
    data LTSURFDIR /1,1,2,2/ 
    integer, dimension(4) :: LTOUT2VEC
    data LTOUT2VEC /-1, +1, +1, -1/
    ! =======================================================================
    
    integer, dimension(2,4) :: LTNEIGH
    data LTNEIGH /-1,0, +1,0, 0,+1, 0,-1/
    integer, dimension(4) :: LTOPST
    data LTOPST /DE, DW, DS, DN/
    
    type BICGLS_MTX
        sequence 
        real(8), pointer, dimension(:,:) :: poiss
        real(8), pointer, dimension(:) :: diag
        integer, pointer, dimension(:,:) :: neigh
        integer :: n, m
    end type

    integer, dimension(4) :: LTS2D ! surf to direction
    integer, dimension(4) :: LTS2N ! surf to normal
    integer, dimension(2,4) :: LTITER
    integer, dimension(2,4) :: LTXS
    integer, dimension(4) :: LTSGF
    integer, dimension(2,4) :: LTMI
    integer, dimension(2,2,4) :: LTFST ! FMC Surface Tally
    integer, dimension(4) :: LTSIGNX ! sign for cross term according to quadrant
    
!    data LTS2D / 2,2, 1,1/
!    data LTS2N / 1,1, 2,2/
!    !data iter /-1,-1, 0,-1, -1,0, 0,0/
!    data LTITER /0,0, -1,0, 0,-1, -1,-1/
!    data LTXS /2,3, 2,4, 1,3, 1,4/
!    data LTSGF /+1, -1, -1, +1/
!    data LTMI /0,0, 1,0, 0,1, 1,1/
!    data LTFST /0,1, 1,1, 0,0, 1,0, 0,0, 0,1, 1,0, 1,1/
!    data LTSIGNX /-1,+1,-1,+1/

    data LTS2D      /1,1, 2,2/
    data LTS2N      /2,2, 1,1/
    data LTITER     /0,-1, -1,-1, 0,0, -1,0/
    data LTXS       /1,3, 1,4, 2,3, 2,4/
    data LTSGF      /-1, +1, +1, -1/
    data LTMI       /0,1, 1,1, 0,0, 1,0/
    data LTFST      /0,0, 0,1, 1,0, 1,1, 0,1, 1,1, 0,0, 1,0/
    data LTSIGNX    /-1,+1, -1,+1/
    integer, parameter :: nring = 10
    
    interface neighf
        module procedure neighfa
        module procedure neighfb
    end interface
contains

subroutine newmtx(mt, n, m)
    type(BICGLS_MTX), intent(inout) :: mt
    integer :: i
    allocate(mt%poiss(m,n))
    allocate(mt%diag(n))
    allocate(mt%neigh(m,n))
    mt%n = n
    mt%m = m
    do i=1, n
        mt%poiss(:,i) = 0
        mt%neigh(:,i) = 0
        mt%diag(i)    = 0
    enddo
end subroutine 
subroutine copymtx(mt2,mt1)
    type(BICGLS_MTX), intent(in) :: mt1
    type(BICGLS_MTX), intent(out) :: mt2
    integer :: i
    if (mt1%m .ne. mt2%m .or. mt1%n .ne. mt2%n) return
    do i=1, mt1%n
        mt2%poiss(:,i) = mt1%poiss(:,i)
        mt2%neigh(:,i) = mt1%neigh(:,i)
        mt2%diag(i)    = mt1%diag(i)
    enddo
end subroutine

function bndryc(oi,oj)
    include "pspec.FI"
    integer :: bndryc, oi, oj
    bndryc = 0
    if (oi<1) then 
        bndryc = DW
    elseif (oj<1) then
        bndryc = DS
    elseif (oi>ncmx) then
        bndryc = DE
    elseif (oj>ncmy) then
        bndryc = DN
    endif
end function 
function bndryf(oi,oj)
    include "pspec.FI"
    integer :: bndryf, oi, oj
    bndryf = 0
    if (oi<1) then 
        bndryf = DW
    elseif (oj<1) then
        bndryf = DS
    elseif (oi>nfmx) then
        bndryf = DE
    elseif (oj>nfmy) then
        bndryf = DN
    endif
end function 

subroutine neighf2(i,j,k,oi,oj)
    include "pspec.FI"
    integer, intent(in) :: i, j, k
    integer, intent(out) :: oi, oj
    oi = i+LTNEIGH(1,k)
    oj = j+LTNEIGH(2,k)
    if (oj .lt. 1 .and. albedo(DN) < 0) oj = nfmy
    if (oj .gt. nfmx .and. albedo(DS) <0) oj = 1
    if (oi .lt. 1 .and. albedo(DW) < 0) oi = nfmx
    if (oi .gt. nfmy .and. albedo(DE) <0) oi = 1
end subroutine

subroutine neighc2(i,j,k,oi,oj)
    include "pspec.FI"
    integer, intent(in) :: i, j, k
    integer, intent(out) :: oi, oj
    oi = i+LTNEIGH(1,k)
    oj = j+LTNEIGH(2,k)
    if (oj .lt. 1 .and. albedo(DN) < 0) oj = ncmy
    if (oj .gt. ncmx .and. albedo(DS) <0) oj = 1
    if (oi .lt. 1 .and. albedo(DW) < 0) oi = ncmx
    if (oi .gt. ncmy .and. albedo(DE) <0) oi = 1
end subroutine

function neighc(i,j,k)
    include "pspec.FI"
    integer, intent(in) :: i,j,k
    integer :: oi, oj 
    integer :: neighc
    neighc=0
    call neighc2(i,j,k,oi,oj)
    if (map_exp2(oi,oj) >0) then
        neighc = i21c(oi,oj)
!        neighc = map_exp(neighc) ! for Vacuum Assembly
    endif
end function 

function neighfa(i,j,k)
    integer, intent(in) :: i,j,k
    integer :: oi, oj 
    integer :: neighfa
    oi = i+LTNEIGH(1,k)
    oj = j+LTNEIGH(2,k)
    if (bndryf(oi,oj)>0) then
        neighfa =0
    else
        neighfa = i21f(oi,oj)
    endif
end function 

function neighfb(i,j,oi,oj,k)
    include "pspec.FI"
    integer, intent(in) :: i,j,k
    integer, intent(out) :: oi, oj 
    integer :: neighfb
    oi = i+LTNEIGH(1,k)
    oj = j+LTNEIGH(2,k)
    neighfb=0
    ! boundary check
    if (bndryf(oi,oj)>0) then
    else
    if (map_exp2(f2cx(oi),f2cy(oj)) > 0) then
        neighfb = i21f(oi,oj)
    else
    endif
    endif
end function 

function i21c(i,j)
    include "pspec.FI"
    integer :: i, j
    integer :: i21c
    i21c = i+(j-1)*ncmx
end function

function i21f(i,j)
    include "pspec.FI"
    integer :: i, j
    integer :: i21f
    i21f = i+(j-1)*nfmx
end function
function inc(i)
    integer :: i, inc 
    i=i+1
    inc = i
end function     


function QuickSearch(array, p, q, k) result(ret)
    real(8), pointer, dimension(:) :: array
    integer :: p, q, k
    integer :: t
    real(8) :: ret
    icall =icall+1
    if (p<=q) then
        t = Partition(Array, p, q)
        if(t .eq. k) then
          ret=Array(t);
          return
        elseif(t > k) then
          ret=QuickSearch(Array, p, t-1, k)
        else
          ret=QuickSearch(Array, t+1, q, k)
        endif    
    endif
contains
function partition(array, p, q) result (ret)

    integer :: ret, p, q
    real(8), pointer, dimension(:) :: array
    
    real(8) :: base
    integer :: i,j
    
    base=array(q)
    i=p-1
    j=p
    do while(j<q)
        if (array(j)<base) then
            i=i+1
            call swap(array, i, j)
        endif
        j=j+1
    enddo 
    i=i+1
    call swap(array, i, q)
    ret=i
    
endfunction

subroutine swap(array, i, j)
    integer :: i, j
    real(8), pointer, dimension(:) :: array
    real(8) :: temp
    temp = array(i)
    array(i)=array(j)
    array(j)=temp
endsubroutine

endfunction

end module
